arXiv:1506.06592vl [astro-ph.EP] 22Jun2015 


Mon. Not. R. Astron. Soc. 000, \T \ -?? (2014) 


Printed 23 June 2015 (MN KTgX style file v2.2) 


Extreme hydrodynamic atmospheric loss near the critical 
thermal escape regime 

N. V. Erkaev, 1,2 , H. Lammer, 3 , P. Odert, 4 , Yu. N. Kulikov, 5 , K. G. Kislyakova 3 

1 Institute for Computational Modelling, Russian Academy of Sciences, Krasnoyarsk 36, Russian Federation 

2 Siberian Federal University, Krasnoyarsk, Russian Federation 

3 Space Research Institute, Austrian Academy of Sciences, Schmiedlstr. 6, A -8Of 2, Graz, Austria 
4 Institute of Physics, University of Graz, Universitdtsplatz 5, A-8010 Graz, Austria 

5 Polar Geophysical Institute, Russian Academy of Sciences, Khalturina Str. 15, Murmansk 183010, Russian Federation 

Released 2014 Xxxxx XX 


ABSTRACT 

By considering martian-like planetary embryos inside the habitable zone of solar-like 
stars we study the behavior of the hydrodynamic atmospheric escape of hydrogen for 
small values of the Jeans escape parameter j3< 3, near the base of the thermosphere, 
that is defined as a ratio of the gravitational and thermal energy. Our study is based 
on a 1-D hydrodynamic upper atmosphere model that calculates the volume heating 
rate in a hydrogen dominated thermosphere due to the absorption of the stellar soft 
X-ray and extreme ultraviolet (XUV) flux. We find that when the /3 value near the 
mesopause/homopause level exceeds a critical value of ~2.5, there exists a steady hy¬ 
drodynamic solution with a smooth transition from subsonic to supersonic flow. For a 
fixed XUV flux, the escape rate of the upper atmosphere is an increasing function of 
the temperature at the lower boundary. Our model results indicate a crucial enhance¬ 
ment of the atmospheric escape rate, when the Jeans escape parameter /3 decreases 
to this critical value. When j3 becomes ^ 2.5, there is no stationary hydrodynamic 
transition from subsonic to supersonic flow. This is the case of a fast non-stationary 
atmospheric expansion that results in extreme thermal atmospheric escape rates. 
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evolution ultraviolet: planetary systems - stars: ultraviolet - hydrodynamics 


1 INTRODUCTION 

Studies related to hydrodynamic escape and evolution of 
planetary atmospheres (e.g., Sekiya et al. 1980a, 1980b, 
1981; Watson et al. 1981; Kasting & Pollack 1983; Yelle 
2004; Tian et al. 2005a, 2005b; Murry-Clay et al. 2009; 
Koskinen et al. 2013; Erkaev et al. 2013, 2014; Lammer et 
al. 2013a, 2014) indicate that the heating of the upper at¬ 
mosphere caused by high XUV radiation and the related 
hydrodynamic expansion of the bulk atmosphere is impor¬ 
tant for the escape of light gases such as hydrogen from 
early planetary atmospheres. The hydrodynamic outflow of 
the atmospheric particles is somewhat similar to that of the 
solar wind described by the well known isothermic model 
of Parker (1964a, 1964b). For escape, the atmospheric gas 
requires a heating source to overcome the gravitational po¬ 
tential of the planetary body. The main heating source in 
planetary atmospheres is provided by the host star in the 
form of XUV radiation and its absorption in the upper at¬ 
mosphere. 

In their pioneering study Watson et al. (1981) proposed 
an analytical hydrodynamic escape model, where the XUV 


heating was assumed to be deposited within a thin layer near 
the planetary surface. In reality the XUV flux is absorbed 
over a wider altitude range (e.g., Tian et al. 2005a; Murry- 
Clay et al. 2009). For instance, Murray-Clay et al. (2009) 
applied a more realistic upper atmosphere heating model 
where the XUV volume heating rate Qxuv can be written 
as 

Qxuv = riI xuv oexp(—r )an, (1) 

with a heating efficiency r/ that corresponds to the fraction of 
the absorbed XUV radiation that is transformed into ther¬ 
mal energy (Shematovich et al. 2014), the XUV absorption 
cross section uxuv which is ~ 5 x 10 -18 cm 2 for a hydrogen 
dominated atmosphere (Beyont & Cairns 1965), the neutral 
atmosphere density n, and optical depth r 



This XUV heating model was further extended by Erkaev 
et al. (2013) by taking into account the angular dependence 
of the heating rate. Generally, the main necessary condition 
for hydrodynamic escape modeling is that the sonic point 
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should be within the collision dominated atmosphere where 
the Knudsen number Kn is sufficiently small. Volkov et al. 
( 2011 ) suggested a criterion where hydrodynamic models 
should be valid without problems as long as the Kn ^0.1. 

In our present study we investigate the escape criteria 
and efficiency of atomic hydrogen from a Mars-like plane¬ 
tary embryo at 1 AU that experiences XUV flux values be¬ 
tween 45-100 times higher than that of today’s solar value. 
A relatively low gravity and temperatures at the base of the 
thermosphere that are ^ 250 K provide low values of the 
Jeans escape parameter /3 


_ GMpim, 
P ~ kTR 


(3) 


which can cause a crucial enhancement of the hydrodynamic 
atmospheric escape rates. With G Newton’s gravitational 
constant, M p \ the mass of a planet, m the mass of the at¬ 
mospheric main species, Boltzmann constant k, atmospheric 
temperature T and planetocentric distance R. For large gas 
giants such as Jupiter the (3 value at the base of the thermo¬ 
sphere is ~ 474 and for hot Jupiter’s such as HD209458b the 
thermosphere base (3 value is ~ 313, while at a distance of 
~ 3-Rpi j3 decreased to ~ 5.6. For low mass planetary bodies 
such as the building blocks of terrestrial planets like martian 
size planetary embryos that outgassed volatiles which pro¬ 
duced steam atmospheres during the solidification of magma 
oceans (Elkins-Tanton 2012), or nebula captured hydrogen 
gas (Lammer et al. 2014), that are exposed by the high XUV 
flux of young host stars inside the habitable zone /? near the 
base of the thermosphere can reach values that are < 3. 

With this in mind, the main aim of our study is to 
investigate the behavior of the hydrodynamic atmospheric 
expansion of a hydrogen dominated upper atmosphere for 
low (3 values. In Sect. 2 we describe the model approach, in 
Sect. 3 we present and discuss our analytical and numerical 
results and their relevance for the evolution of initial water 
inventories of terrestrial planets. 


2 MODEL EQUATIONS 

2.1 Hydrodynamic upper atmosphere model 


To study the XUV-heated upper atmosphere structure and 
thermal escape rates of exposed hydrogen atoms, we apply a 
1-D hydrodynamic upper atmosphere model with an energy 
absorption model, that is described in detail in Erkaev et al. 
(2013, 2014) and Lammer et al. (2013a, 2014). The model 
solves the system of the hydrodynamic equations for mass, 


dpR 2 dpvR 2 
dt + dR 


(4) 


momentum, 


dpvR 2 9 [ R 2 {pv 2 + P)] 

dt OR 

and energy conservation, 


pgR 2 + 2 PR, 


(5) 


dR 2 


pv 2 I p 

2 ^ ( 7 -l) 


dt. 


+ 


dvR 2 


pv 2 I ~jP 
2 + ( 7 - 1 ) 


OR 


pvR 2 g + QxuvP 2 - 


( 6 ) 


The distance R corresponds to the radial distance from the 
centre of the planet, p, P, T, v are the mass density, pressure, 


temperature and velocity of the XUV exposed and heated 
non-hydrostatic outward flowing bulk atmosphere. 7 is the 
polytropic index, g the gravitational acceleration and Qxuv 
is the XUV volume heating rate (Erkaev et al. 2013, 2014) 


ir/2+a 


s(l/r) 


Qxuv = »70-xuvn/xuv; 


J(r,6)xuvsm(6)d0, (7) 


where r is the normalized radial distance r = R/Ro, Ro is 
the lower boundary radius at the base of the thermosphere 
near the mesopause/homopause level, 7 is the heating ef¬ 
ficiency which is typically 15% (Shematovich et al. 2014), 
7xuv is the XUV flux at the upper boundary, crxuv is the 
XUV absorption cross section in a hydrogen atmosphere 
(Chassefiere 1996; Erkaev et al. 2013; Lammer et al. 2014), 
and J(r,6)x uv is a dimensionless function describing the 
variation the XUV flux with respect to the radial distance 
due to the atmospheric absorption, 


J(r,0)x uv = exp 




r 2 sin(<9))~ 1 / 2 £d£ 


( 8 ) 


Here R* is the upper boundary radius, n = n/no, where 
no is the normalized density at the lower boundary, and 
a = axvvnoRo is a dimensionless constant parameter. For¬ 
mula © describes the XUV flux intensity as a function of 
spherical coordinates r and 9. In the particular case when 
9=0, this formula can be simplified to the expression 

J(r, 9)x uv = exp(-r), (9) 

where r is given by equation ©. Therefore, our heating 
model is be quite similar to that of Murray-Clay et al. (2009) 
for the central radial direction with zero spherical angle. 

I 11 case of a stationary outflow regime, we can neglect 
the time derivatives and obtain the ordinary system of equa¬ 
tions as follows 


pv 




J_ 

2 ^ ( 7 _i)p 


dv ( 
dR + « 

+ $ 


dR 


pvR . 2 = vk, 

( 10 ) 


( 11 ) 

7~ ~ p dR’ 

= QxvvR 2 , 

( 12 ) 


where $ is the gravitational potential, and 4/ is the escaping 
mass flux per unit steradian. 

Integrating equation m we get a relationship between 
the escape rate and the total absorbed energy. 


jf QxvvR 2 dR 

[v* 2 /2 + AJ? — vo 2 /2 — 7/(7 — l)(fcTo/m)]' 


(13) 


Here Vo and v* are the flow velocities at the lower and upper 
boundaries. The total absorbed energy (energy deposition) 
is related to the incoming XUV flux 


[R, 

47t J QxvvR 2 dR = r/Ixuv^ReS, 


(14) 


where R e g is the effective XUV absorbtion radius 


R e fi = Ro 


1+ / (1- J(s,0))2sds 


1/2 


(15) 
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Introducing definition £ = mA<5>/(fcTo), and using (1141) we 
rewrite equation m 


A = 


4i r4/ = 
1 


ylxUV'fReS^ 

A4> 


A, 


[l + (n* 2 -no 2 )/(2A$)- 7 /[( 7 - 1)£]] ’ 


(16) 

(17) 


For large values of the Jeans escape parameter £, the effec¬ 
tive radius R e s is close to the base of the thermosphere Ro , 
the ratio n* 2 /A4? can then be neglected and the coefficient 
A m 1. In such a case equation © yields the well known 
energy limited escape formula 


47r4' = 


^XUV^eg 

AT 


(18) 


In cases of a small £ value near Ro, the escape rate be¬ 
comes much larger than that predicted by the energy lim¬ 
ited escape formula in because the effective radius (J? e ff) 
increases substantially, and in equation m the terms kT/m 
and AT become comparable to each other. These terms be¬ 
come equal to each for a critical value of £ 


mA$/(ff 0 ) = /3 C = 


7-1 


(19) 


In such a case the enthalpy of the lower atmospheric gas is 
large enough to provide expansion of the atmosphere against 
gravity, and an additional heating is not needed for that. The 
escape rate is expected to have a pile-up when the Jeans 
escape parameter /3 approaches to /3 C . 


3 ANALYTICAL AND NUMERICAL RESULTS 


The system of equations (11011111121) can be transformed to 
a pair ordinary differential equations with respect to the 
normalized entropy and velocity 


dV 

dr 


V 

r 


dS_ _ (7 - 1 )qr 2( ' y - 1) 
dr V(2 —r)r(7-i) ’ 

'/ 3 / r - 2 7 S (^)' y " 1 + ( 7 - l )^ 

: 7 s(^r l -H ' 


( 20 ) 

( 21 ) 


Here T is the normalized escape rate 
4/ 

r =- . -, (22) 

(no^kTo/mRl) 

S is the entropy-like function, S = P /n r_1 , P is the normal¬ 
ized pressure, P = P/(nokTo)), V is the normalized velocity 
v = v/y/m/ m, r is the normalized distance, r = R/Ro, 
and q is the normalized heating rate per one particle 


QxxiVy/rnRo 

q ~ n(fcTo) 3 / 2 ' ( j 

The denominator of equation (ED means a difference be¬ 
tween the sonic speed squared and the bulk velocity squared. 
At the sonic point, where the denominator vanishes, the nu¬ 
merator must vanish as well. In the subsonic region the de¬ 
nominator becomes positive, and thus the numerator should 
be also positive in the most part of the region below the 
sonic point in order to provide an acceleration of the flow. 
However, near the lower boundary the numerator can be 


negative in case of a low /?. Therefore, for a low /3 the nu¬ 
merator has to change its sign in the subsonic region and 
at the sonic point. In such a case, the velocity behavior is 
not monotonic. First it decreases near the lower boundary, 
and then it starts to increase until the sonic point. There 
are two ways for finding the escape rate and radial profiles 
of the velocity, density and temperature. 

The first way is to apply the so-called shooting method 
for solving the ordinary differential equations with the 
boundary condition at the sonic point. This method will 
allow us to determine the unknown escape rate T. For the 
particular unique value of the escape rate it is possible to 
pass by the sonic point from the subsonic lower boundary 
to the supersonic upper boundary. 

The second way is to use the time relaxation method 
by solving the non-steady system of the hydrodynamic equa¬ 
tions. We find the numerical solution by a finite difference 
numerical scheme described in the previous publications of 
Erkaev et al. (2013). Steady-state profiles are obtained by 
time relaxation of the numerical solution. We used both 
methods mentioned above and obtained similar results. In 
particular we investigate the behavior of the solution when 
the escape parameter £ becomes close to the critical value 
7/(7 - 1 ). 

The results shown in Figs. 1-3 and Table 1 are based 
on an assumed hydrogen dominated upper atmosphere that 
originated either from a hot steam atmosphere or nebula 
captured hydrogen gas above a planetary body with the size 
and mass of Mars in an orbit inside the habitable zone of 
a Sun-like star. The high XUV fluxes of young solar-like 
stars will most likely dissociate hydrogen bearing molecules 
such as H 2 O or H 2 and in the upper atmosphere so that it 
should be mainly dominated by atomic hydrogen (Kastin- 
gand & Pollack 1983; Chassefiere 1996; Yelle 2004; Lammer 
2013). During such protoatmosphere conditions one can ex¬ 
pect that the surface is much hotter (i.e. solidified magma 
ocean, magma ocean in mush stage, surface bombarded fre¬ 
quently by planetesimals) compared to that of present Mars 
in the Solar System (e.g., Marcq 2013; Lebrun et al. 2013; 
Erkaev et al. 2014; Lammer et al. 2014). Because of this, 
the results of 1-D radiativeconvective atmospheric models 
that studied the coupling between hot planetary surfaces 
and surrounding steam atmospheres (Marcq 2012) indicate 
that the mesopause/homopause level zq of low mass bodies 
such as Mars may extend from ~100 km as on present day 
Mars up to distances that can be located a few 100 Kilome¬ 
ter higher. For this reason we assume our lower boundary 
conditions at a distance Ro = R p 1 + zo with zo = 433 km, 
with a mesopause/homopause temperature that correspond 
to the equilibrium or skin temperature at 1 AU, assumed to 
be To=247 K (£=2.9), 256 K (£=2.8), and 266 K (/3=2.7). 
For the number density at that level we assume a similar 
value no = 5 x 10 12 cm -3 (e.g., Tian et al. 2005b; Erkaev et 
al. 2013, 2014; Lammer et al. 2013, 2014). 

On the left hand side, Fig. 1 shows the analytical Mach 
number profiles corresponding to three different values of 
velocity at the lower boundary and a Jeans escape param¬ 
eter £ of 2.7. On the right hand side, this figure shows the 
analytical velocity profiles, and also the numerical profile 
for £=2.7. The numerical one is obtained by time relaxation 
method based on the finite difference numerical scheme de¬ 
scribed in Erkaev et al. (2013). The analytical curves are 
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Figure 1. Left hand side: Analytical sonic Mach number profiles for /3 =2.7. Right hand side: Analytical and numerical velocity profiles 
for /3 =2.7; here curves (a),(b),and (c) correspond to different initial velocity values: 0,03843, 0.03841 and 0.0383, respectively. 



Figure 2. Analytical Mach number profiles (left) and velocity 

obtained by the integration of equations (HU EU) for three 
different initial velocity conditions: Vo = 0.03843, 0.03841 
and 0.0383. The true initial value Vo = 0, 03843 corresponds 
to a separatrix, which provides acceleration of particles from 
subsonic to supersonic values. This separatrix is rather close 
to the profile obtained by the numerical relaxation method. 

Fig. 2 shows the analytical sonic Mach number and ve¬ 
locity profiles corresponding to a lower j3 parameter value of 
2.5. In this case (/3=2.5), there is no acceleration at all, when 
the lower boundary velocity is subsonic. By considering dif¬ 
ferent subsonic velocity conditions, we obtain only decreas¬ 
ing velocity functions. For any subsonic velocity conditions 
at the lower boundary, the velocity does not have monotonic 
increase until supersonic values. Our results show that the 
escape rate increases enormously, when the escape parame¬ 
ter /? approaches its critical value. 

Fig. 3 shows the normalized escape rate (per sterad.) 
and the total escape rate (per sec.) as functions of the Jeans 
escape parameter /3 for XUV flux values that are 45 and 
100 times higher compared to that of today’s Sun. For a 
quite small decrease of the /3 parameter, the escape rate in¬ 
creases by a factor of 10. However, the sonic point position, 
and the corresponding velocity and temperature have a mi¬ 
nor change when the escape rate is approaching the critical 



;s (right) corresponding to different initial velocities for /3 = 2.5 

value. An exception is the density at the sonic point, which 
has a substantial increase proportionally to the escape rate. 

Table 1 presents the total escape rate, effective radius, 
and also hydrodynamic parameters such as the flow veloc¬ 
ity, density and temperature corresponding to two partic¬ 
ular points. The first one ( R s ) is the sonic point where 
the sonic Mach number is ~1. The second one (R*) is the 
point where the Knudsen number is ~0.1. The model re¬ 
sults are based of the Jeans escape parameter near the 
mesopause/hompause level for three values near the critical 
one, namely 2.7, 2.8, and 2.9 and two XUV flux intensities 
that are 45 and 100 times higher than that of today’s solar 
value. 

One can see that with our assumed input parameters, 
hydrodynamic conditions are valid far beyond R s is reached. 
The escape rate rises up to ~3 times if the /3 value near 
the critical one, changes only slightly from 2.9 to 2.7. Thus, 
extreme uncontrolled blow-off will occur for smaller values. 

Chamberlain (1963) and Watson et al. (1981) described 
/3 = 1.5 as a critical value for atmospheric “blow-off” corre¬ 
sponding to the case when the thermal energy is equal to the 
gravitational energy. In the model of Parker (1964a, 1964b), 
the critical (5 value is equal to 2. For lower values of /?, a hy¬ 
drodynamic flow cannot have transition from subsonic to a 


© 2014 RAS, MNRAS 000, [[]-?? 


































5 




Figure 3. Normalized escape rate (per sterad.) and the total escape rate (per sec.) as functions of the Jeans escape parameter 0 for 
XUV flux values that are 45 and 100 times higher compared to that of today’s Sun. 


Table 1. Main parameters characterizing the hydrodynamic escape flow from a martian-like body at 1 AU, surrounded by a hydrogen 
envelope (i.e., dissociated H2O molecules, captured hydrogen from a planetary nebula: R s is the sonic point distance, and R* is the 
distance where the Knudsen number is ~0,1. The 0 values of 2.9, 2.8 and 2.7 correspond to To values neat the mesopause/homopause 
level of 247 K, 256 K, and 266 K, respectively. 


XUV 

Escape rate, 10 33 s 1 

Rea/Ro 

Ra/Ro 

V 0 /V To 

Na/N 0 

Ta/T 0 

R*/Ro 

V*/V To 

N* /No , 10" 9 

T* /To 

0 

45 

128 

13 

40 

0,31 

70 

0,06 

8000 

2,6 

0,29 

1,4 

2,7 

45 

64 

10 

40 

0,31 

35 

0,06 

4500 

2,1 

0,48 

0,9 

2,8 

45 

40 

9 

39 

0,32 

34 

0,06 

3000 

1,9 

0,66 

0,7 

2,9 

100 

199 

13 

31 

0,37 

143 

0,08 

8000 

3,5 

0,25 

2,3 

2,7 

100 

105 

10 

31 

0,37 

81 

0,08 

5000 

2,9 

0,46 

1,7 

2,8 

100 

63 

9 

30 

0,38 

54 

0,08 

3300 

2,5 

0,66 

1,3 

2,9 


supersonic regime. A continuous acceleration is possible only 
in the supersonic regime. This means that the flow velocity 
should be supersonic already at the lower boundary. 

From our hydrodynamic model we find the critical /3 
parameter is ~2.5, which corresponds to a condition that 
the enthalpy is equal to the gravitational energy. It is also 
worth to mention that our results agree with the kinetic 
simulations of Volkov et al. (2011) whose simulations yield 
a very sharp growth of the escape rate, when the 0 param¬ 
eter decreases from 2.7 to 2. For 0 values ^ 2 this kinetic 
simulations predict only supersonic flow. For larger 0 values 
the kinetic model of Volkov et al. (2011) predicts only sub¬ 
sonic flow because it does not have heating above Ro because 
the incoming XUV radiation was assumed to be absorbed 
around Ro- Or model includes a strong XUV heating above 
Ro, and therefore we have acceleration from a subsonic to a 
supersonic regime for /3 > 2.5. 

Small 0 values can be reached on smaller planetary bod¬ 
ies such as planetary embryos after magma oceans solidified 
(e.g.,Lammer et al. 2013b; Lebrun et al. 2013), or lower mass 
hydrogen-dominated exoplanets in close orbital distances. 
One can see from Table 1, Mars-type bodies with hydro¬ 
gen dominated upper atmospheres will reach this critical 
conditions quite easy. Smaller bodies will reach the criti¬ 
cal hydrodynamic escape regime even easier. Under extreme 
conditions the atmosphere at the base of the thermosphere 
may be hot so that a supersonic flow can occur around the 
mesopause/homopause level. Such a very fast non-steady ex¬ 
pansion of an atmosphere may act like an explosion that en¬ 
hances the loss of water from the building blocks of the ter¬ 
restrial planets, that may finally accrete drier as expected. 


H 2 O and other volatiles may be stepwise outgassed and 
potentially lost to a great extent, as illustrated in Fig. 10 in 
Lammer et al. (2013b), before accretion ends. If wet plane¬ 
tary embryos lose much of their initial H 2 O by fast hydro- 
dynamic escape during their growth to the final planetary 
body, the volatile content that is outgassed in the final stage 
would be lower than expected. Such a scenario would agree 
with the hypothesis of Albarede and Blichert- Toft (2007), 
which is that the terrestrial planets in the Solar System ac¬ 
creted dry and obtained most of their water during the late 
veneer via impacts. 

However, the fast expansion will be connected by strong 
adiabatic cooling which has a feedback on the 0 parameter, 
that will grow above the critical value where a stationary 
hydrodynamic escape regime will most likely be established. 


4 CONCLUSION 

Analyzing the hydrodynamic atmospheric escape solutions 
for a hot martian-type planetary body (i.e., large planetary 
embryos in an inner planetary system, etc.), we found a 
critical value for the Jeans escape parameter 0 C ~ 2.5 cor¬ 
responding to a bifurcation of the atmospheric expansion 
regime. A hydrodynamic flow acceleration from subsonic to 
supersonic velocities is possible for (3 values that exceed this 
critical value. The true value of the escape rate corresponds 
to the separatrix velocity profile, which goes from the lower 
boundary through the sonic point towards infinity. The es¬ 
cape rate is found to increase crucially, when 0 approaches 
the critical value. When /3 is below the critical value, the 
velocity becomes supersonic even at the lower boundary. 
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Regarding this, the following scenario is expected. A su¬ 
personic flow at the lower thermosphere means a very fast 
non-steady expansion of the atmosphere like an explosion. 
This expansion might cause a strong adiabatic cooling of the 
atmosphere. Such process might result in a decrease of the 
temperature, and consequently an increase of the /3 param¬ 
eter. Finally the /3 parameter will grow above the critical 
value, and thus a stationary hydrodynamic escape can be 
established. 
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